Multifractal Scaling in the Bak-Tang-Wiesenfeld Sandpile and Edge Events 
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An analysis of moments and spectra shows that, while 
the distribution of avalanche areas obeys finite size scaling, 
that of toppling numbers is universally characterized by a 
full, nonlinear multifractal spectrum. Rare, large avalanches 
dissipating at the border influence the statistics very sensibly. 
Only once they are excluded from the sample, the conditional 
toppling distribution for given area simplifies enough to show 
also a well defined, multifractal scaling. The resulting picture 
brings to light unsuspected, novel physics in the model. 

PACS numbers: 64.60. Lx, 64.60.Ak, 05.40. +j, 05.65.+b 

Finite size scaling (FSS) 0] is a widely adopted frame- 
work for the description of finite, large systems near crit- 
icality. In the last decade, after the work of Bak et al. 
Q , much attention has been devoted to a class of models 
in which criticality is spontaneously generated by the dy- 
namics itself, without the necessity of tuning parameters. 
This self-organized criticality (SOC) has been advocated 
as a paradigm for a wide range of phenomena, from earth- 
quakes to interface depinning, economics and biological 
evolution The prototype model of SOC is the two 
dimensional (2D) Bak, Tang, and Wiesenfeld sandpile 
(BTW) , which represents a system driven by a slow 
external influx, dissipated at the borders through a local, 
nonlinear mechanism. 

In spite of its apparent simplicity and relative ana- 
lytical tractability g-Q], the 2D BTW resisted, so far, 
all theoretical attempts to fully and exactly characterize 
its scaling . These attempts were essentially all based 
on the FSS ansatz. Numerical approaches, also based on 
FSS |8| , led to rather scattered and sometimes contradic- 
tory numerical results |)] [ll| , which hardly concile with 
existing theoretical conjectures. Thus, with its intriguing 
intractability, BTW scaling remains a formidable chal- 
lenge for nonequilibrium statistical mechanics j[2| and it 
is very important to check if FSS works in this context. 

In this Letter we apply a new strategy of data collec- 
tion and interpretation, in order to determine to what 
extent the FSS ansatz can be applied, or rather has to 
be modified, for a correct description of the 2D BTW. 
Our results are striking and largely unexpected: while 
compelling evidence is obtained that the probability dis- 
tribution functions (pdf) of some quantities obey FSS, for 
other magnitudes, whose fractal dimensions can widely 
fluctuate within the the nonlinear dynamics, this is def- 
initely not the case. Following our protocol of analy- 



sis, we demonstrate that the well known difficulties in 
the description of the BTW are due to unexpected, very 
nonstandard features of its dynamical behavior. In the 
BTW, relations between different key quantities do not 
reduce to standard power laws, as in FSS, and are sub- 
stantially influenced by the infrared cutoff given by the 
size of the system. The peculiar fluctuations characteriz- 
ing intermittent dissipation at the cutoff scale, provide a 
dynamical mechanism for unusual deviations from finite 
size, and even multifractal scaling. 

We consider the BTW on a square lattice box A; to 
each site i we associate an integer height Z\ > 0, the 
number of grains. When z% exceeds a threshold z c = 4, 
site i topples: z, — > Zj — 4, while for the nearest neighbors 
j of i, Zj — > Zj + 1. At the boundary less than 4 neigh- 
bor sites are upgraded with consequent grain dissipation. 
Further instabilities can be created by the first toppling. 
An avalanche is the set of the s topplings necessary to 
reach a stable system configuration after addition of one 
grain ( Zk — > Zk + 1 at some randomly chosen k s A), a 
is the number of lattice sites toppling at least once dur- 
ing the avalanche. A sequence of avalanches is created 
by successive random additions. After sufficiently many 
grains, thanks to dissipation at the borders, the sandpile 
reaches a steady state. We analyzed up to 10 8 avalanches 
in this state for L = 128, 256, 512 and 1024. 

The pdf for a and s do not reveal characteristic scales 
intermediate between lattice spacing and linear pile size 
L. FSS postulates for them the forms 

P.(s,L) = 8- T 'F B (a/L D ') (1) 
P a (a ) L) = a- T *F a (a/L D *) ) 

and usually assumes that different quantities characteriz- 
ing an avalanche are simply related by power laws, deter- 
mined by sharply peaked conditional pdf. For example, 
one expects that, given a, s ~ a 7 with D s — jD a , 
and D s (D a ) representing the fractal dimension of s (a). 
Without sticking to the FSS form (|l|), scaling of P s is 
most generally described by the multifractal spectrum 

Ha) _MZ%*m. (2) 

log(L) 

/ is the Legendre transform of the moment scaling func- 
tion a (q) defined by: 
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(s q ) L = I P s {s,L) S ods~L^\ (3) 

i.e. a (q) = sup a [qa + /(«)]. Analogous definitions apply 
to the spectrum g(f3) (/3 = log a/ log L) and the moment 
exponent p(q) of P a . If Eqs. ([!]) hold, / (a) = — (r s — \)a 
for < a < D s and / = — oo for a > D s . Consistently 
ct (q) = D s (q — t s + 1) for q > t s — 1 and a (q) = for 
q < t s — X. Corresponding expressions for et-quantities 
hold if FSS is valid. So, within FSS, both / and ct are 
piece- wise linear functions of their arguments [ p"3"| . 

A reliable way to establish if Eqs.([y) hold, is by check- 
ing the above linearity of ct and p in significant ranges 
of q. Being extrapolated for L — ► oo from finite L 
moments, p and ct provide a very asymptotic charac- 
terization of the pdf's. While for P a a constant gap 
&p(q) = P(q + 1) - p(o) ~ 2.02 ± 0.03 - D a establishes 
already for q = 1 (Table §), for P s , ct'(1) ~ 2.5, while Act 
steadily increases from Act(1) ~ 2.70 to Act(8) ~ 2.92. 
Act(1) ~ 2.7, was also found in Ref. p5[, based on large L 
data. Thus, unlike for P a , and in violation of FSS, for P s 
there is clearly no constant gap in the range q > I. The 
gap tends to rise to ~ 3.0 for increasing q. Therefore, we 
should also expect f(a) > — oo as long as a < 3.0. 

Fig.|l| reports g and / as obtained by the data collapse 
technique in Ref. ||. The linear form g (/3) = — (r a — 1)(3 
is well verified for (3 <■ 1.5 with an estimated slope 
—0.19 ± 0.01. In the region (3 > 1.5, the collapse gets 
worse. One expects D a = 2 and, in fact, the poor g- 
collapse for (3 < 2 is consistent with the infinite disconti- 
nuity of a FSS spectrum with D a — 2: curves for various 
L smooth out to different degrees such discontinuity, and 
underestimate g for (3 < 2. Assuming a linear g in the 
whole domain < (3 < 2 , and r a = 1/5 exactly, as sug- 
gested by the estimated initial slope, {o)l should scale 
with p(l) = supp [g{(3) + [3] = -2/5 + 2 = 1.6, in nice 
agreement with our determination p(l) = 1.59 ± 0.02. 
Evidence of FSS for P a comes also from the fact that 
a standard FSS collapse (r Q = 6/5; D a — 2) works 
very well (Fig. ([!]), inset). For a < 2 the collapsed / 
is very close to linear and overlaps with the expected g 

(/(2) 0.39). Thus, an acceptable FSS form of P s 

should assume r s = t q , in order to be consistent with 
the well collapsed, initial part of the plots. Within such 
assumption, D s = 2.5 would be imposed by the exact 
result ct(1) = 2 [§ (we find ct(1) = 1.99 ± 0.02). In- 
deed, sup Q [/(a) + a] = 2 in such case, the sup being 
attined at a = D s = 2.5. The hypothetical linear spec- 
trum should therefore have support in < a < 2.5, and 
satisfy / (5/2) = -1/2. In Ref. @, a relatively limited 
analysis of ct suggested Act(<7) ~ 2.5 for q > 1. Such con- 
stant gap would leave room to a FSS approximation for 
P s , with r s ~ 6/5 and D s ~ 2.5, and a linear ct deviating 
from the measured one possibly only at very low g's. Ac- 
cording to Eq. (0), for such effective P s one would also 
have D s (2 — t s ) ~ 2 = ct(1). However, this does not agree 
with the plots which show that / > — oo at least up to 
a ~ 3.0. Indeed, curves for various L collapse rather well 



for a < 3.0, and clearly suggest a support of / asymptot- 
ically bounded by a ~ 3.0, consistent with the trend of 
Act (Table |). This bound follows also from the leftward 
trend of the curves for increasing L in the region a ^ 3, 
where collapse gets worse. An even approximate t s ex- 
ponent, so extensively discussed in the last decade, can 
not be simultaneously consistent with the initial slope 
1/5 of /, Act(oo) - 3 and ct(1) = 2. Full consis- 
tency, in particular with ct(1) = 2, can be recovered by 
assuming that / is indeed linear, with the same slope 
— 1/5 as g, up to a = 2.5 (/(2.5) = —1/2), but has a 
nonlinear continuous drop in the range 2.5 < a < 3.0. 
The slight underestimation by the plot for a ~ 2.5 ( 
/(5/2) ~ —0.57) should again be imputed to roundoff 
and slower L-convergence in correspondence to the ma- 
jor bending of / [ p7| . 

The striking difference described above between P s and 
P a suggests that the conditional pdf, C, such that 

P s (s,L)= J daC(s\a,L)P a (a,L), (4) 

should have unusual structure. Within FSS one would 
assume C ~ 6(s — a 1 ) for a < L 2 and L — > oo. 
Such C would lead to 7 = (t s — l)/(r a — 1). Thus, 
t s = r a would imply 7=1, while D s = 2.5 would give 
7 = D s /D a = 1.25, already contradicting FSS. In fact 
C is a complex, broad pdf. Its properties elucidate and 
confirm our conclusions for P a and P s . Fig. |^ reports val- 
ues of a = logs/ log L vs. (3 = log a/ log L for L = 1024 
avalanches. Ratios 7 = a/(3 range between 1 ( s > a) 
and ~ 1.25, and their spread is not modified appreciably 
by sampling data for progressively larger L's. For f3 < 2 
also some a/ (3 > 1.25 are found, which are too rare to 
be displayed in Fig. ||. If a relation s ~ a 1 would hold, 
data should coalesce into a straight line with slope 7. To 
the contary, points are quite spread and form an open 
angle, rather than a narrow strip, as is instead the case, 
e.g., for the radius of gyration, r, of the surface covered 
by the avalanche (Fig. ||) . 

C-moments can be assumed to scale as (s q ) a = 
f ds C(s\a,L)s q - a K ^'?) = . One would hope 

exponents like k to be well defined, i.e. independent of 
(3, for L — > 00. Furthermore, with a FSS C oc 8{s — a 7 ), 
one should find (3k{(3, q) j q — [3j independent of q. In 
Fig.^we plot ^((s 9 )^ 9 )/ log L versus (3 for various q and 
L = 1024. The curves, which correspond to (3K,(P,q)/q 
asymptotically, do not overlap. Moreover, for each q the 
plots have pronounced curvature, especially for (3 > 1.7. 
This curvature does not decrease appreciably, or increase 
signalling crossovers, if one considers progressively larger 
L's. Thus, in spite of the sensible spread of the various 
curves in Figj3| (which also persists for increasing L), the 
complex scaling of C-moments can not even be classified 
as multifractal. Indeed, in that case one should still re- 
quire (3k j q to be linear in (3 |lj]. On the other hand, the 
/3-nonlincarity described above is essential in view of our 
conclusions on g. Since (s 9 )l = J da P a (a, L) (s q ) a , we 
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must have a(q) = sup g[g ((3) +f3n(q, 0)}. With g = — (3/5, 
as argued above, the sup for q — 1 should fall at /3 = 2 
(Fig. ||), and be equal to — 2/5+2k(1, 2). Thus, one must 
find K[l, 2) = 1.2. If k(1, (3) would remain constant with 
its initial value 1.04 ± 0.02 [[l8| in the whole /3-range, we 
could certainly not find <r(l) = 2. Our determinations 
consistently extrapolate to 2) ~ 1.20 (Fig. ||). 

The /3-dependence of re is crucial for the consistency 
of the overall scaling properties and is determined by 
large dissipating avalanches. These are intermittent, edge 
events, which guarantee grain conservation at stationar- 
ity pq] . Remarkably, if only nondissipative avalanches 
are sampled, one obtains a different conditional pdf, Co, 
whose re is now independent of (3 (Fig.||). Thus, for Co 
it makes sense to discuss a multifractal spectrum 



(3h(j) 



lp g(/^ dsC (s\a,L)) 
log(L) 



(5) 



As a function of 7, h measures the density of points 
{a = P^f, (3) along a narrow vertical strip at fixed f3 in 
the plot correspoding to that in Fig.||. Fig.|| reports a 
data collapse of ^.-curves at /3 = 1.6 [^9|. Even if the 
collapse suffers of poor L-asymptoticity, all curves cross 
nicely at a point corresponding to 7 ~ 1.25, where their 
trend for increasing L inverts itself. Collapses at different 
/3's give very similar results. Thus, 7 ~ 1.25 qualifies as a 
maximum exponent for nondissipating, bulk avalanches. 
An independent, accurate determination based on rank 
ordering statistics gave 7 = 1.28 ± 0.03. 

7 = 5/4 was conjectured as a unique exponent for all 
avalanches in Ref. M. Different values Jig] , sometimes 
larger than 1.25 |15|7 were conjectured in the recent lit- 
erature on the basis of FSS and numerical results. Figs. 
I and I show that the determination of 7 is not even a 
well defined task, and can be discussed, within a multi- 
fractal framework, only once edge avalanches are elimi- 
nated. Indeed, for nondissipative bulk avalanches, 7 has 
a /3-independent, broad, nonlinear spectrum, h, with 1 as 
most probable, and ~ 1.25 as extreme, most rare values. 

The above analysis identifies the inadequacies of the 
scaling theory in Ref. ||. Crucial to that approach was 
the assumption of a relation, n oc a' 27-2 )/ 2 = a 1 / 4 , sat- 
isfied by n, the maximum number of topplings in an 
avalanche (realized at the injection point), and a, seen 
as the area of the first in a succession of waves || . We 
found that also the 71-pdf at given a is broad, and 1/4 
is only the approximate maximum value of the exponent 
realized by bulk sandpile dynamics. 

The constant gap Aer ~ 2.5 at high q postulated in 
Ref. |l6| would have been compatible with a narrow, 
or even point-like spectrum for dissipating avalanches 
(a ~ 2.5). To the contrary, it turns out that also dis- 
sipating avalanches are multifractal (2.0 < a < 3.0) con- 
sistent with our more correct estimate of Act at high q. 
Moreover, a frequency ~ with £ = | is precisely 
associated to dissipating avalanches with a > 2.50. 

In summary, we showed how the long standing puzzle 



of 2D BTW scaling finds a solution in a strong viola- 
tion of FSS by both P s and C. While P a obeys FSS, 
with D a = 2 and r a = 6/5 as most plausible exponents, 
nonlinear multifractal spectra are needed to character- 
ize P s and Co- Our results throw light on the intrigu- 
ing difficulty of this model. In spite of the several ex- 
actly known steady state properties, the belief that the 
2D BTW should be "easily" solvable appears unjusti- 
fied. The unusual scaling pattern discovered here, which 
is not found in more simple systems, like directed sand- 
piles p2| , constitutes genuine novel physics and enhances 
the paradigmatic role of the BTW for statistical mechan- 
ics out of equilibrium. 

The most striking feature of this pattern is the bend- 
ing of the curves in FigJ|, determined by the effect of 
intermittent, edge avalanches, and essential in order to 
fulfill the exact constraint of local, Laplacian conserva- 
tion at stationarity (cr(l) = 2). In fact the very presence 
of edge avalanches allows the existence of a peculiar, bulk 
multifractal scaling of Co- The first moment of the s-pdf 
of nondissipative avalanches scales ~ L 18 , rather than 
~ L 2 , as conservation imposes in the case of P s . The role 
played here by intermittent edge avalanches has striking 
analogies with the intermittent phenomena in fully devel- 
oped turbulence, where they cause deviation from pure 
Kolmogorov scaling (^o|. Results of the present analysis 
have been used most recently to elucidate quantitative 
connections between BTW scaling and the notoriously 
difficult problem of turbulence in 3D j2lj| . 
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q 


1 


2 


3 


4 


5 


6 


7 


Ap 


2.02 


2.02 


2.02 


2.03 


2.02 


2.02 


2.02 


Aa 


2.70 


2.83 


2.87 


2.91 


2.91 


2.92 


2.92 



TABLE I. Gaps at different q for P a and P s . The estimated 
uncertainty is ±0.03. 




Iog(s)/log(0.3 L) log(a)/log(L) 

FIG. 1. Multifractal collapses for g and /; a standard FSS 
collapse of P > (a,L) is in the inset. The subscript > indicates 
integrated pdf. 
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(3=log(a)/log(L) 

FIG. 2. (/3, a) points for avalanches with L — 1024 (upper); 
similar plot for (/3, log r / log L) (lower). The slope in the latter 
case is D" 1 = 0.50 ±0.01. 
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FIG. 3. The lower set of curves refers to all avalanches, 
while the upper one (shifted by 1 along y axes) pertains to 
the nondissipative ones, q — 1, 4, 10, 16 moving upwards. 
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